threshold <- function(x, condition = 0, comparator = "<") {

        if (comparator == ">") {
            ifelse(x > condition, return(1), return(0))
        } else if (comparator == ">=") {
            ifelse(x >= condition, return(1), return(0))
        } else if (comparator == "<=") {
            ifelse(x <= condition, return(1), return(0))
        } else if (comparator == "<") {
            ifelse(x < condition, return(1), return(0))
        } else if (comparator == "==") {
            ifelse(x == condition, return(1), return(0))
        } else if (comparator == "!=") {
            ifelse(x != condition, return(1), return(0))
        } else {
            return(FALSE)
        }

}

threshold_less_than_five <- function(x) {
    return(threshold(x, condition = 5, comparator = "<"))
}

threshold_less_than_zero <- function(x) {
    return(threshold(x, condition = 0, comparator = "<"))
}

threshold_greater_than_zero <- function(x) {
    return(threshold(x, condition = 0, comparator = ">"))
}

threshold_greater_than_one <- function(x) {
    return(threshold(x, condition = 1, comparator = ">"))
}

is_NA <- function(x) {
    ifelse(is.na(x) || is.nan(x), return(0), return(x))
}
# Question 1
palo <- readr::read_csv("data/uscities.csv") %>%
    st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
    filter(city == "Palo", state_name == "Iowa") %>%
    st_transform(5070)

palo_bbox <- palo %>%
    st_buffer(5000) %>%
    st_bbox() %>%
    st_as_sfc() %>%
    st_as_sf()

# Question 2.2 (src: R/scenes.R)
# To build the palo-flood-scene.csv, uncomment below:
# source("R/scenes.R")

# palo_scene <- readr::read_csv("data/palo-flood-scene.csv")[1,]$download_url %>%
#     lsat_scene_files()

# palo_scene_urls <- palo_scene %>%
#     filter(grepl(paste0('B', 1:6, ".TIF$", collapse = "|"), palo_scene$file)) %>%
#     arrange(file) %>%
#     pull(file)

# Question 2.3

# lsat_image() was not working for me due to a partition error/encryption
# in my HDD, which gave me the error "Invalid cross-device link" when
# file.rename() was called. So, I worked around this by manually
# downloading the *.TIF files via `wget` into the ~/.cache/landsat-pds/...
# directory (In Fedora).
# st <- sapply(palo_scene_urls, lsat_image)

# This results in *.TIF files still being recognized in the cache:
st <- lsat_cache_list()

bands <- stack(st) %>%
    setNames(paste0("band", 1:6))

bands_crop <- bands %>%
    crop(st_transform(palo_bbox,"+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs"))

bands_crop <- setNames(bands_crop, c("Coastal", "Blue", "Green", "Red", "NIR", "SWIR 1"))

Landsat Band Attributes

The attributes of our stacked bands are given as,

  • nrow(): 7,811
  • ncol(): 7,681
  • ncell(): 59,996,291
  • nlayers(): 6
  • CRS: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
  • Resolution: 30 x 30

The attributes of our cropped bands are given as,

  • nrow(): 340
  • ncol(): 346
  • ncell(): 117,640
  • nlayers(): 6
  • CRS: +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs
  • Resolution: 30 x 30

Plotting Bands & Applying Stretches

RGB (4,3,2)

Traditional Color Infrared (5,4,3)

False Color, Water Focus (5,6,4)

False Color, Agriculture Focus (6,5,2)

Why Stretch?

Stretching in this case is a way of increasing clarity by normalizing the distribution of brightness across a raster. In general, these diagrams from NeonScience display the general concept behind stretching:

Darker Lighter

When stretching, there are two methods: Linear or Histogram. Linear is the method shown in the diagrams above, where points are taken and linearly scaled. Histogram stretching occurs conceptually by stretching the ends of the intesity (brightness/contrast) of a raster: Histogram

NVDI <- (bands_crop[[5, ]] - bands_crop[[4, ]]) / (bands_crop[[5, ]] + bands_crop[[4, ]])
NVDI_mask <- calc(NVDI, threshold_less_than_zero) %>%
    calc(is_NA)

NDWI <- (bands_crop[[3, ]] - bands_crop[[5, ]]) / (bands_crop[[3, ]] + bands_crop[[5, ]])   
NDWI_mask <- calc(NDWI, threshold_greater_than_zero) %>%
    calc(is_NA)

MNDWI <- (bands_crop[[3, ]] - bands_crop[[6, ]]) / (bands_crop[[3, ]] + bands_crop[[6, ]])
MNDWI_mask <- calc(MNDWI, threshold_greater_than_zero) %>%
    calc(is_NA)

WRI <- (bands_crop[[3, ]] + bands_crop[[4, ]]) / (bands_crop[[5, ]] + bands_crop[[6, ]])
WRI_mask <- calc(WRI, threshold_greater_than_one) %>%
    calc(is_NA)

SWI <- (1 / (bands_crop[[2, ]] - bands_crop[[6, ]])^0.5) %>%
    calc(is_NA)
SWI_mask <- calc(SWI, threshold_less_than_five) %>%
    calc(is_NA)

water_features <- stack(c(NVDI, NDWI, MNDWI, WRI, SWI)) %>%
    setNames(c("NVDI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(water_features, col = colorRampPalette(c("blue", "white", "red"))(256))

water_features_mask <- stack(c(NVDI_mask, NDWI_mask, MNDWI_mask, WRI_mask, SWI_mask)) %>%
    setNames(c("NVDI", "NDWI", "MNDWI", "WRI", "SWI"))
plot(water_features_mask, col = colorRampPalette(c("white", "blue"))(256))

set.seed(123456)
kmeans_floods <- getValues(water_features_mask) %>%
    na.omit() %>%
    kmeans(10)

kmeans_raster <- water_features$NVDI
values(kmeans_raster) <- kmeans_floods$cluster